Introduction

In this document, we will outline the Bayesian analogs of the statistical analyses described in lecture 1 (Github code).

Load packages

library(dplyr)
library(tibble)
library(purrr)
library(tidyr)
library(forcats)
library(gtools)
library(patchwork)
library(broom)
library(broom.mixed)
library(modelr)
library(brms)
library(tidybayes)
library(ggdist)
library(bayesplot)
library(ggplot2)
library(knitr)
BRM_BACKEND <- ifelse(require("cmdstanr"), 'cmdstanr', 'rstan')

Dataset

The following dataset is from experiment 2 of “How Relevant are Incidental Power Poses for HCI?” (Jansen & Hornbæk, 2018). Study participants were asked to either make an expansive posture or a constrictive posture before performing a task. The experiment investigated whether posture could potentially have an effect on risk taking behavior.

First, we load the data.

pose_df = readr::read_csv("data/poses_data.csv", show_col_types = FALSE) %>%
  mutate(condition) %>%
  group_by(participant)

The data has been aggregated for each participant: - condition = expansive indicates expansive posture, and condition = constrictive indicates constrictive posture - The dependent variable is change which indicates the percentage change in risk-taking behavior. Thus, it is a continuous variable.

For the purposes of this demo, we are only concerned with these two variables. We can ignore the other variables for now.

head(pose_df %>% select(participant, condition, change))
participant condition change
1 expansive 3.362832
2 constrictive 29.147982
3 expansive 25.409836
4 constrictive 54.069767
5 expansive -36.644592
6 constrictive 29.756098

Intuition of Bayesian Statistics

The Bayesian t-test (BEST) assumes that the data in the two conditions arises from two separate t-distributions. In the following section, we will describe the process for one of the conditions in the experiment.

We will use the \(Normal(\mu = 20, \sigma = 20)\) as the prior distribution.

First, we define some functions for manual calculation of the posterior normal distribution:

sigma_post = function(sigma_prior, sigma, n = 1) {
  sqrt(1 / (1 / (sigma_prior^2) + n / (sigma^2)))
}
mu_post = function(mu_prior, sigma_prior, mu, sigma, n = 1) {
  tau = sigma_post(sigma_prior, sigma, n)
  (tau^2 / sigma_prior^2)*mu_prior + (n * tau^2 / sigma^2)*mu
}
d.p2 = tibble(
  group = c("prior", "expansive", "constrictive", "posterior"), 
  mu = c(20, 32.82, 31.61, mu_post(20, 20, 32.82, 7.52)), 
  sd = c(20, 7.52, 7.06, sigma_post(20, 7.52))
) %>%
  mutate(
    cutoff_group = list(c(1:7)),
    cutoff = list(c(0, 15, 25, 30, 32.82, 40, 100))
  ) %>%
  unnest(c(cutoff_group, cutoff)) %>%
  mutate(
    cutoff = if_else(group == 'prior', 100, cutoff),
    x = map(cutoff, ~ seq(from = -40, ., by = 0.1)),
    y = pmap(list(x, mu, sd), ~ dnorm(..1, ..2, ..3))
  )

In the plot below, we show the raw data distribution for the two conditions:

p1 = pose_df %>%
  ggplot() +
  geom_point(aes(x = change, y = condition, colour = condition), 
             position = position_jitter(height = 0.1), alpha = 0.7) +
  scale_color_theme() + 
  labs(y = "Condition") +
  theme(
    legend.position = "none", 
    axis.line.y = element_blank(), 
    axis.ticks.y = element_blank(),
    axis.title.y = element_blank()
  ) +
  scale_x_continuous(limits = c(-150, 250), breaks = seq(-150, 250, by = 50))

p2.blank = tibble(y = c("expansive", "constrictive"), x = 0) %>%
  ggplot(aes(x, y)) +
  scale_color_theme() + 
  theme_density

cowplot::plot_grid(p2.blank, p1, nrow = 2)

First, we define a function to help us plot different cutoff groups.

plot_preliminary <- function(data, group_num, group_name){
  data %>%
  filter(cutoff_group == group_num & group %in% group_name) %>%
  unnest(c(x, y)) %>%
  ggplot(aes(x, y)) +
  #geom_line(aes(color = group), size = 1) +
  # density
  geom_area(aes(fill = group, color = group), position = "identity", linewidth = 1, alpha = 0.3) +
  scale_x_continuous(limits = c(-40, 100)) +
  scale_y_continuous(limits = c(0, 0.1)) +
  scale_color_theme() +
  theme_density
}
cowplot::plot_grid(
  plot_preliminary(d.p2, 0, NULL), 
  p1, nrow = 2)

Next, we plot the prior density:

cowplot::plot_grid(
  plot_preliminary(d.p2, 1, 'prior'), 
  p1, nrow = 2)

Then we describe step by step, how the likelihood is computed:

cowplot::plot_grid(
  plot_preliminary(d.p2, 1, 'expansive'), 
  p1, nrow = 2)

cowplot::plot_grid(
  plot_preliminary(d.p2, 2, 'expansive'), 
  p1, nrow = 2)

cowplot::plot_grid(
  plot_preliminary(d.p2, 3, 'expansive'), 
  p1, nrow = 2)

cowplot::plot_grid(
  plot_preliminary(d.p2, 5, 'expansive'), 
  p1, nrow = 2)

cowplot::plot_grid(
  plot_preliminary(d.p2, 6, 'expansive'), 
  p1, nrow = 2)

cowplot::plot_grid(
  plot_preliminary(d.p2, 7, 'expansive'), 
  p1, nrow = 2)

We want to compute the posterior, which is the product of the prior and likelihood:

cowplot::plot_grid(
  plot_preliminary(d.p2, 7, c('prior', 'expansive', 'posterior')), 
  p1, nrow = 2)

Of course, if you have gganimate and magick working on your computer, you can try the following code to get a gif! Here, we are not going to show the code.

Fumeng: the content below needs better documentation….

Model 1: Equal standard deviations

Step 1: model specification

Fumeng: I hope to keep this example simple. Is there a simple way to justify the choice of \(\nu\)?

Here we use a mathematical expression for the model above. Fumeng: need explanations here…

\[ \begin{align} y_{i} &\sim \mathrm{Student\_t}(\mu, \sigma_{0}, \nu_{0}) \\ \mu &= \beta_{0} + \beta_{1} * x_i \\ \sigma_{0} &\sim \mathrm{HalfNormal}(0, 10) \\ \beta_{0} &\sim \mathrm{?} \\ \beta_{1} &\sim \mathrm{Normal}(0,2) \\ \nu_{0} &\sim \mathrm{?} \\ i & \in \{\mathrm{expansive}, \mathrm{constrictive}\} \end{align} \]

We translate this thought into brms formula using the function bf.

model.1.formula <- bf(
                      # we think change is affected by different conditions.
                      change ~ condition,
                      # to tell brms which response distribution to use
                      family = student()
                    )

Now we can use get_prior() to inspect the available priors and formula. As what we learned, we have priors for

tibble(get_prior(model.1.formula, pose_df))
prior class coef group resp dpar nlpar lb ub source
b default
b conditionexpansive default
student_t(3, 22.6, 38.8) Intercept default
gamma(2, 0.1) nu 1 default
student_t(3, 0, 38.8) sigma 0 default

prior checks

Look at the default priors

cowplot::plot_grid(
tibble(x = qstudent_t(ppoints(n = 500), df = 3, mu = 22.6, sigma = 38.8)) %>% 
  ggplot() +
  geom_density(aes(x = x), fill = theme_yellow, color = NA) +
  ggtitle('student_t(3, 22.6, 38.8) for Intercept') +
  coord_cartesian(xlim = c(-300, 300),expand = c(0)),
tibble(x = qgamma(ppoints(n = 500), shape = 2, rate = .1)) %>% 
  ggplot() +
  geom_density(aes(x = x), fill = theme_yellow, color = NA) +
  ggtitle('Gamma(2, 0.1) for nu') +
  coord_cartesian(xlim = c(1, 100), expand = 0),
tibble(x = qstudent_t(ppoints(n = 500), df = 3, mu = 0, sigma = 38.8)) %>% 
  ggplot() +
  geom_density(aes(x = x), fill = theme_yellow, color = NA) +
  ggtitle('student_t(3, 0, 38.8) for sigma') +
  coord_cartesian(xlim = c(0, 400),expand = c(0)),
ncol = 3)

prior check for default priors

model.1.checks_default <- brm(
  model.1.formula,
  data = pose_df, 
  family = student_t(),
  prior = c(
    prior(normal(0, 3), class = 'b')
  ),
  sample_prior = 'only',
  #backend = BRM_BACKEND,
  file = 'rds/model.1.checks_default.rds',
  file_refit = 'on_change' 
)

plot

model.1.defaultpriorsamples <- 
  model.1.checks_default %>% 
    epred_draws(tibble(condition = c('expansive', 'constrictive')))
head(model.1.defaultpriorsamples)
condition .row .chain .iteration .draw .epred
expansive 1 NA NA 1 19.425813
expansive 1 NA NA 2 38.452771
expansive 1 NA NA 3 -23.169166
expansive 1 NA NA 4 48.733025
expansive 1 NA NA 5 72.826955
expansive 1 NA NA 6 -2.607648 
model.1.defaultpriorsamples %>% 
  ggplot(aes(x = .epred, group = condition)) +
  geom_density(alpha = .5, color = NA, adjust = 2, fill = theme_yellow) +
  theme_density_x + 
  ggtitle('Checks default priors of model.1')

Our priors

cowplot::plot_grid(
tibble(x = qnorm(ppoints(n = 1000), mean = 0, sd = 2)) %>% 
  ggplot() +
  geom_density(aes(x = x), fill = theme_yellow, color = NA) +
  ggtitle('Normal(0, 2) for b & Intercept') +
  coord_cartesian(xlim = c(-8, 8), expand = 0),
tibble(x = qnorm(ppoints(n = 1000), mean = 0, sd = 10)) %>% 
  ggplot() +
  geom_density(aes(x = x), fill = theme_yellow, color = NA) +
  ggtitle('HalfNormal(0, 10) for sigma') +
  coord_cartesian(xlim = c(0, 50), expand = c(0)),
ncol = 2)

model.1.checks <- brm(
  model.1.formula,
  data = pose_df, 
  family = student_t(),
  prior = c(
    prior(normal(0, 2), class = 'b'),
    prior(normal(0, 10), class = 'sigma', lb = 0)
  ),
  sample_prior = 'only',
  #backend = BRM_BACKEND,
  file = 'rds/model.1.checks.rds',
  file_refit = 'on_change' 
)

plot

model.1.priorsamples <- 
  model.1.checks %>% 
    epred_draws(tibble(condition = c('expansive', 'constrictive')))

Compare two sets of priors

cowplot::plot_grid(
  
model.1.defaultpriorsamples %>% 
  ggplot(aes(x = .epred, group = condition)) +
  geom_density(alpha = .5, color = NA, adjust = 2, fill = theme_yellow) +
  theme_density_x + 
  #coord_cartesian(xlim = c(-800, 800)) + 
  ggtitle('Checks default priors of model.1')
,
model.1.priorsamples %>% 
  ggplot(aes(x = .epred, group = condition)) +
  geom_density(alpha = .5, color = NA, adjust = 2, fill = theme_yellow) +
  theme_density_x + 
  #coord_cartesian(xlim = c(-800, 800)) + 
  ggtitle('Checks our priors of model.1')
,
ncol = 1)

Step 2: model fitting

model.1 <- brm(
  model.1.formula,
  data = pose_df, 
  family = student_t(),
  prior = c(
    prior(normal(0, 2), class = 'b'),
    prior(normal(0, 10), class = 'sigma', lb = 0)
  ),
  backend = BRM_BACKEND,
  file = 'rds/model.1.rds'
)

Step 3: check posteriors

aspect 1: mcmc traces

color_scheme_set("teal")
mcmc_trace(model.1, facet_args = list(ncol = 4))

aspect 2: model metrics

Rhat, ESS, etc here

summary(model.1)
##  Family: student 
##   Links: mu = identity; sigma = identity; nu = identity 
## Formula: change ~ condition 
##    Data: pose_df (Number of observations: 80) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept             22.84      4.49    14.50    31.77 1.00     2686     2092
## conditionexpansive    -0.26      1.97    -4.10     3.62 1.00     2726     2987
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    28.35      4.13    20.76    36.92 1.00     2180     2520
## nu        3.68      2.35     1.57     9.04 1.00     2172     2158
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

aspect 3: visually

Fumeng: do we need this?

model.1.predictions <- 
  model.1 %>% 
    predicted_draws(tibble(condition = c('expansive', 'constrictive')))
head(model.1.predictions)
condition .row .chain .iteration .draw .prediction
expansive 1 NA NA 1 33.21591
expansive 1 NA NA 2 19.94745
expansive 1 NA NA 3 47.89480
expansive 1 NA NA 4 19.25483
expansive 1 NA NA 5 40.19085
expansive 1 NA NA 6 -14.57711 
plot_predictions <- function(model, df = NULL, title = ''){
  
  if(is.null(df))
      df = tibble(condition = c('expansive', 'constrictive'), 
                           participant = c(-1, -1))
  model %>% 
    predicted_draws(df,
                    seed = 1234,
                    ndraws = NULL,
                    allow_new_levels = TRUE,
                    sample_new_levels = 'uncertainty') %>% 
    ggplot(aes(x = .prediction, fill = condition)) +
    geom_density(alpha = .5, size = 1, adjust = 2, color = NA) +
    scale_color_theme() +
    theme_density_x + 
    scale_y_continuous(breaks = 0, labels = 'constrictive') + 
    scale_x_continuous(breaks = seq(-150, 250, by = 50)) +
    coord_cartesian(xlim = c(-150, 250)) + 
    ggtitle(paste0('Posterior predictions of ', deparse(substitute(model))))
}
cowplot::plot_grid(plot_predictions(model.1), 
                   p1 + 
                    scale_x_continuous(breaks = seq(-150, 250, by = 100)) +
                    coord_cartesian(xlim = c(-150, 250)), 
                   nrow = 2)

model.1.posteriors <- 
  model.1 %>% 
    epred_draws(tibble(condition = c('expansive', 'constrictive')))
plot_posteriors <- function(model, df = NULL, title = ''){
  
  if(is.null(df))
     df = tibble(condition = c('expansive', 'constrictive'))
  
  model %>% 
    epred_draws(df, 
                # ignoring random effects if there is any
                #seed = 1234,
                ndraws = NULL,
                re_formula = NA) %>% 
    ggplot(aes(x = .epred, fill = condition)) +
    geom_density(alpha = .5, size = 1, adjust = 2, color = NA) +
    scale_color_theme() +
    theme_density_x + 
    scale_y_continuous(breaks = 0, labels = 'constrictive') + 
    scale_x_continuous(limits = c(-150, 250), breaks = seq(-150, 250, by = 50)) +
    ggtitle(paste0('Posterior means of ', deparse(substitute(model))))
}
cowplot::plot_grid(plot_posteriors(model.1), p1, nrow = 2)

Model 2: the BEST test model

Step 1: model specification

This model is the BEST test model as described by Kruschke in the paper Bayesian estimation supersedes the t-test. In this model, \(\beta\) indicates the mean difference in the outcome variable between the two groups (in this case, the percent change in the BART scores). We fit different priors on \(\beta\) and set different weights on these priors to obtain our posterior estimate.

\[ \begin{align} y_{i} &\sim \mathrm{T}(\nu, \mu, \sigma) \\ \mu &= \beta_{0} + \beta_{1} * x_i \\ \sigma &= \sigma_{a} + \sigma_{b}*x_i \\ \beta_{1} &\sim \mathrm{Normal}(\mu_{0}, \sigma_{0}) \\ \sigma_a, \sigma_b &\sim \mathrm{Cauchy}(0, 2) \\ \nu &\sim \mathrm{exp}(1/30)\\ i & \in \{\mathrm{expansive}, \mathrm{constrictive}\} \end{align} \]

model.2.formula <- bf(# we think change is affected by different conditions.
                      change ~ condition,
                      sigma ~ condition,
                      # to tell brms which response distribution to use
                      family = student())
tibble(get_prior(model.2.formula, pose_df))
prior class coef group resp dpar nlpar lb ub source
b default
b conditionexpansive default
student_t(3, 22.6, 38.8) Intercept default
gamma(2, 0.1) nu 1 default
b sigma default
b conditionexpansive sigma default
student_t(3, 0, 2.5) Intercept sigma default

prior check

Fumeng: not sure qcauchy has the correct paramaterization

cowplot::plot_grid(
  tibble(x = qnorm(ppoints(n = 1000),  mean = 0, sd = 2)) %>% 
  ggplot() +
  geom_density(aes(x = x), fill = theme_yellow, color = NA) +
  ggtitle('Normal(0 ,2) for Intercept') +
  coord_cartesian(expand = c(0)),
tibble(x = qexp(ppoints(n = 1000), rate = 0.0333)) %>% 
  ggplot() +
  geom_density(aes(x = x), fill = theme_yellow, color = NA) +
  ggtitle('exponential(0.0333) for nu') +
  coord_cartesian(expand = 0),
tibble(x = qcauchy(ppoints(n = 1000), location = 0, scale = 2)) %>% 
  ggplot() +
  geom_density(aes(x = x), fill = theme_yellow, color = NA) +
  ggtitle('cauchy(0, 2) for sigma') +
  coord_cartesian(expand = c(0)),
ncol = 3)

model.2.priorchecks <- brm(
  model.2.formula,
  data = pose_df, 
  family = student_t(),
  prior = c(
    prior(normal(0, 2), class = 'b'),
    prior(cauchy(0, 2), class = 'b', dpar = 'sigma'),
    prior(exponential(0.0333), class = 'nu')
  ),
  sample_prior = 'only',
  #backend = BRM_BACKEND,
  file = 'rds/model.2.priorchecks.rds',
  file_refit = 'on_change' 
)
model.2.priorchecks %>% 
 epred_draws(tibble(condition = c('expansive', 'constrictive'))) %>% 
  ggplot(aes(x = .epred, group = condition)) +
  geom_density(alpha = .5, color = NA, adjust = 2, fill = theme_yellow) +
  theme_density_x + 
  ggtitle('Checks priors of model.2')

Step 2: model fiting

model.2 <- brm(
  model.2.formula,
  data = pose_df, 
  family = student_t(),
  prior = c(
    prior(normal(0, 2), class = 'b'),
    prior(cauchy(0, 2), class = 'b', dpar = 'sigma'),
    prior(exponential(0.0333), class = 'nu')
  ),
  backend = BRM_BACKEND,
  file = 'rds/model.2.rds',
  file_refit = 'on_change' 
)
## Running MCMC with 4 sequential chains...
## 
## Chain 1 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 1 Iteration:  100 / 2000 [  5%]  (Warmup) 
## Chain 1 Iteration:  200 / 2000 [ 10%]  (Warmup) 
## Chain 1 Iteration:  300 / 2000 [ 15%]  (Warmup) 
## Chain 1 Iteration:  400 / 2000 [ 20%]  (Warmup) 
## Chain 1 Iteration:  500 / 2000 [ 25%]  (Warmup) 
## Chain 1 Iteration:  600 / 2000 [ 30%]  (Warmup) 
## Chain 1 Iteration:  700 / 2000 [ 35%]  (Warmup) 
## Chain 1 Iteration:  800 / 2000 [ 40%]  (Warmup) 
## Chain 1 Iteration:  900 / 2000 [ 45%]  (Warmup) 
## Chain 1 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 1 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
## Chain 1 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
## Chain 1 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
## Chain 1 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
## Chain 1 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
## Chain 1 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
## Chain 1 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
## Chain 1 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
## Chain 1 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 2 Iteration:  100 / 2000 [  5%]  (Warmup) 
## Chain 2 Iteration:  200 / 2000 [ 10%]  (Warmup) 
## Chain 2 Iteration:  300 / 2000 [ 15%]  (Warmup) 
## Chain 2 Iteration:  400 / 2000 [ 20%]  (Warmup) 
## Chain 2 Iteration:  500 / 2000 [ 25%]  (Warmup) 
## Chain 2 Iteration:  600 / 2000 [ 30%]  (Warmup) 
## Chain 2 Iteration:  700 / 2000 [ 35%]  (Warmup) 
## Chain 2 Iteration:  800 / 2000 [ 40%]  (Warmup) 
## Chain 2 Iteration:  900 / 2000 [ 45%]  (Warmup) 
## Chain 2 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 2 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
## Chain 2 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
## Chain 2 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
## Chain 2 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
## Chain 2 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
## Chain 2 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
## Chain 2 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
## Chain 2 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
## Chain 2 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
## Chain 2 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 2 finished in 0.1 seconds.
## Chain 3 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 3 Iteration:  100 / 2000 [  5%]  (Warmup) 
## Chain 3 Iteration:  200 / 2000 [ 10%]  (Warmup) 
## Chain 3 Iteration:  300 / 2000 [ 15%]  (Warmup) 
## Chain 3 Iteration:  400 / 2000 [ 20%]  (Warmup) 
## Chain 3 Iteration:  500 / 2000 [ 25%]  (Warmup) 
## Chain 3 Iteration:  600 / 2000 [ 30%]  (Warmup) 
## Chain 3 Iteration:  700 / 2000 [ 35%]  (Warmup) 
## Chain 3 Iteration:  800 / 2000 [ 40%]  (Warmup) 
## Chain 3 Iteration:  900 / 2000 [ 45%]  (Warmup) 
## Chain 3 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 3 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 3 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
## Chain 3 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
## Chain 3 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
## Chain 3 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
## Chain 3 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
## Chain 3 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
## Chain 3 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
## Chain 3 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
## Chain 3 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
## Chain 3 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 3 finished in 0.1 seconds.
## Chain 4 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 4 Iteration:  100 / 2000 [  5%]  (Warmup) 
## Chain 4 Iteration:  200 / 2000 [ 10%]  (Warmup) 
## Chain 4 Iteration:  300 / 2000 [ 15%]  (Warmup) 
## Chain 4 Iteration:  400 / 2000 [ 20%]  (Warmup) 
## Chain 4 Iteration:  500 / 2000 [ 25%]  (Warmup) 
## Chain 4 Iteration:  600 / 2000 [ 30%]  (Warmup) 
## Chain 4 Iteration:  700 / 2000 [ 35%]  (Warmup) 
## Chain 4 Iteration:  800 / 2000 [ 40%]  (Warmup) 
## Chain 4 Iteration:  900 / 2000 [ 45%]  (Warmup) 
## Chain 4 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 4 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 4 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
## Chain 4 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
## Chain 4 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
## Chain 4 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
## Chain 4 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
## Chain 4 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
## Chain 4 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
## Chain 4 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
## Chain 4 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
## Chain 4 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.9 seconds.

Step 3: posterior checks

summary(model.2)
##  Family: student 
##   Links: mu = identity; sigma = log; nu = identity 
## Formula: change ~ condition 
##          sigma ~ condition
##    Data: pose_df (Number of observations: 80) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                   24.48      4.69    15.61    34.16 1.00     2536
## sigma_Intercept              3.41      0.20     3.00     3.80 1.00     2430
## conditionexpansive          -0.14      1.97    -3.94     3.67 1.00     4222
## sigma_conditionexpansive     0.15      0.22    -0.29     0.57 1.00     3731
##                          Tail_ESS
## Intercept                    2052
## sigma_Intercept              2331
## conditionexpansive           2869
## sigma_conditionexpansive     2860
## 
## Family Specific Parameters: 
##    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## nu     6.16      9.42     1.73    24.18 1.00     1848     1380
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
wrap_plots(plot_predictions(model.2) + scale_x_continuous(limits = c(-200,250), expand = c(0,0), breaks = seq(-150, 250, by = 50)) 
                   , 
                   p1 + scale_x_continuous(limits = c(-200,250), expand = c(0,0), breaks = seq(-150, 250, by = 50)) , 
                   nrow = 2)

# draws.model.3 = 
wrap_plots(
  plot_posteriors(model.2), 
  p1,
nrow = 2)

model.2.posteriors <- model.2 %>% 
                epred_draws(tibble(condition = c('expansive', 'constrictive')), 
                #seed = 1234,
                ndraws = NULL,
                re_formula = NA)

Model 3: the BEST model with better intercepts

Step 1: model specification

\[ \begin{align} y_{i} &\sim \mathrm{T}(\nu, \mu, \sigma) \\ \mu &= \beta_{i,j} + \beta_{1} * x_i \\ \sigma &= \sigma_{a} + \sigma_{b}*x_i \\ \beta_{1} &\sim \mathrm{Normal}(0, 2) \\ \sigma_a, \sigma_b &\sim \mathrm{HalfNormal}(0, 10) \\ \nu &\sim \mathrm{exp}(1/30)\\ i & \in \{\mathrm{expansive}, \mathrm{constrictive}\}\\ j & \in \{1, ..., \mathrm{N}\} \end{align} \]

model.3.formula <- bf(# we think change is affected by different conditions.
                      change ~ condition + (1|participant),
                      sigma ~ condition,
                      # to tell brms which response distribution to use
                      family = student())
tibble(get_prior(model.3.formula, pose_df))
prior class coef group resp dpar nlpar lb ub source
b default
b conditionexpansive default
student_t(3, 22.6, 38.8) Intercept default
gamma(2, 0.1) nu 1 default
student_t(3, 0, 38.8) sd 0 default
sd participant default
sd Intercept participant default
b sigma default
b conditionexpansive sigma default
student_t(3, 0, 2.5) Intercept sigma default

Step 2: model fiting

model.3 <- brm(
  model.3.formula,
  data = pose_df, 
  family = student_t(),
  prior = c(
    prior(normal(0, 2), class = 'b'),
    prior(normal(0, 2), class = 'b', dpar = 'sigma'),
    prior(normal(0, 2), class = 'sd', group = 'participant', lb = 0),
    prior(exponential(0.0333), class = 'nu')
  ),
  #backend = BRM_BACKEND,
  file = 'rds/model.3.rds',
  file_refit = 'on_change' 
)
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include   -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## 
## SAMPLING FOR MODEL '098ecfbfafd8b9d1ae94bd6b9b0d7ae1' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 3.6e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.36 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.196147 seconds (Warm-up)
## Chain 1:                0.180298 seconds (Sampling)
## Chain 1:                0.376445 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '098ecfbfafd8b9d1ae94bd6b9b0d7ae1' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 1.4e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.1924 seconds (Warm-up)
## Chain 2:                0.183196 seconds (Sampling)
## Chain 2:                0.375596 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '098ecfbfafd8b9d1ae94bd6b9b0d7ae1' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 1.8e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.18 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.222635 seconds (Warm-up)
## Chain 3:                0.182754 seconds (Sampling)
## Chain 3:                0.405389 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '098ecfbfafd8b9d1ae94bd6b9b0d7ae1' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 1.4e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.209288 seconds (Warm-up)
## Chain 4:                0.179886 seconds (Sampling)
## Chain 4:                0.389174 seconds (Total)
## Chain 4:

Step 3: posterior checks

summary(model.3)
##  Family: student 
##   Links: mu = identity; sigma = log; nu = identity 
## Formula: change ~ condition + (1 | participant) 
##          sigma ~ condition
##    Data: pose_df (Number of observations: 80) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Group-Level Effects: 
## ~participant (Number of levels: 80) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.59      1.22     0.04     4.53 1.00     2896     1951
## 
## Population-Level Effects: 
##                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                   24.56      5.03    14.96    34.83 1.00     5232
## sigma_Intercept              3.41      0.20     3.00     3.79 1.00     5833
## conditionexpansive          -0.09      2.00    -4.06     3.84 1.00    10152
## sigma_conditionexpansive     0.15      0.22    -0.29     0.57 1.00     8612
##                          Tail_ESS
## Intercept                    3332
## sigma_Intercept              3450
## conditionexpansive           2853
## sigma_conditionexpansive     2961
## 
## Family Specific Parameters: 
##    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## nu     5.98      7.52     1.70    25.11 1.00     5049     2838
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
cowplot::plot_grid(plot_predictions(model.3), p1, nrow = 2)

cowplot::plot_grid(plot_posteriors(model.3), p1, nrow = 2)

model.3.posteriors <-
  model.3 %>% 
  epred_draws(tibble(condition = c('expansive', 'constrictive')),
               re_formula = NA)
wrap_plots(
model.3.posteriors %>%
  mutate(model = 'model 3') %>%
  rbind(
    model.2.posteriors %>% 
      mutate(model = 'model 2'))  %>% 
  rbind(
    model.1.posteriors %>% 
      mutate(model = 'model 1')) %>% 
ggplot() + 
  geom_density(aes(x = .epred, fill = condition), adjust = 1.5, color = NA, alpha = .5) +
  facet_grid(model ~ .) +
  scale_x_continuous(limits = c(-150, 250), breaks = seq(-150, 250, by = 50)) +
  scale_color_theme() +
  theme_density_x +
  ggtitle('Compare the means of all three models'),
p1, nrow = 2, heights = c(4,1.5))

Model 4: Poisson / Negative-Binomial Regression model

Step 1: model specification

\[ \begin{align} y_{i} &\sim \mathrm{Neg-Binomial}(\mu, \phi) \\ \mu &= \beta_{i,j} + \beta_{1} * x_i \\ \beta_{1} &\sim \mathrm{Normal}(0, 1) \\ i & \in \{\mathrm{expansive}, \mathrm{constrictive}\}\\ j & \in \{1, ..., \mathrm{N}\} \end{align} \]

pose_raw_df = read.csv("data/posture_data-raw.csv") %>%
  mutate(participant = factor(participant)) %>%
  rename(trial = trial.number)

head(pose_raw_df)
participant condition total.money trial trial.money exploded pumps life
1 expansive 910 0 62 0 62 72
1 expansive 910 1 0 1 27 27
1 expansive 910 2 47 0 47 98
1 expansive 910 3 0 1 86 86
1 expansive 910 4 60 0 60 104
1 expansive 910 5 0 1 26 26 
model.4 = brm(pumps ~ 1 + condition + trial + (trial | participant),
                      data = pose_raw_df, family = negbinomial,
                      prior = c(prior(normal(0, 10), class = Intercept),
                                prior(normal(0, 1), class = b),
                                prior(gamma(0.01, 0.01), class = shape)
                      ),  # this is the brms default
                      file = 'rds/model.4.rds',
                      file_refit = 'on_change' ,
                      iter = 4000, warmup = 1000, cores = 4, chains = 4)
prior_summary(model.4)
prior class coef group resp dpar nlpar lb ub source
normal(0, 1) b user
b conditionexpansive default
b trial default
normal(0, 10) Intercept user
lkj_corr_cholesky(1) L default
L participant default
student_t(3, 0, 2.5) sd 0 default
sd participant default
sd Intercept participant default
sd trial participant default
gamma(0.01, 0.01) shape 0 user 
summary(model.4)
##  Family: negbinomial 
##   Links: mu = log; shape = identity 
## Formula: pumps ~ 1 + condition + trial + (trial | participant) 
##    Data: pose_raw_df (Number of observations: 2400) 
##   Draws: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
##          total post-warmup draws = 12000
## 
## Group-Level Effects: 
## ~participant (Number of levels: 80) 
##                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)            0.38      0.04     0.31     0.45 1.00     4207
## sd(trial)                0.01      0.00     0.00     0.01 1.00     3694
## cor(Intercept,trial)    -0.62      0.12    -0.82    -0.35 1.00     9370
##                      Tail_ESS
## sd(Intercept)            6443
## sd(trial)                4265
## cor(Intercept,trial)     7921
## 
## Population-Level Effects: 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept              3.52      0.06     3.41     3.63 1.00     2993     5184
## conditionexpansive    -0.04      0.07    -0.18     0.10 1.00     3431     5201
## trial                  0.01      0.00     0.00     0.01 1.00     8485     8646
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape     7.67      0.28     7.15     8.22 1.00    17805     8817
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Step 2: model interpretation

The results of this model are on the log-odds scale. What do the coefficients mean? The simplest way is to simply transform the data into a more interpretable scale and visualise the results:

# draws.model.4 = 
wrap_plots(
  plot_predictions(model.4, df = crossing(condition = c('expansive', 'constrictive'), 
                                          trial = 0:29,
                                          participant = 0)), 
  p1,
nrow = 2)

We generate the average of 30 trials.

model.4.posteriors <-  model.4 %>% 
    epred_draws(crossing(condition = c('expansive', 'constrictive'),     
                       trial = 0:29,
                       participant = 0),
                allow_new_levels  = TRUE) %>% 
  mutate(model = 'model 4')  %>% 
  group_by(.draw, condition) %>% 
  summarise(.epred = mean(.epred))
wrap_plots(
  model.4.posteriors %>% 
 ggplot(aes(x = .epred, fill = condition)) +
    geom_histogram(alpha = .5, size = 1, binwidth = 5, color = NA) +
    scale_color_theme() +
    facet_wrap(condition ~ ., ncol = 1) + 
    theme_density_x + 
    scale_x_continuous(limits = c(-100, 250), breaks = seq(-100, 250, by = 50)) +
    ggtitle(paste0('Posterior means of ', deparse(substitute(model)))), 
  p1,
nrow = 2)

From this, we can see that there does not appear to be a difference between the two conditions.

Reporting

Let’s use model 4

model

We will need the ESSs and RHat from this print. Also, if you want to report the standard deviation of random intercepts or CIs for other paramaters. You can find them via summary(..)

summary(model.4)
##  Family: negbinomial 
##   Links: mu = log; shape = identity 
## Formula: pumps ~ 1 + condition + trial + (trial | participant) 
##    Data: pose_raw_df (Number of observations: 2400) 
##   Draws: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
##          total post-warmup draws = 12000
## 
## Group-Level Effects: 
## ~participant (Number of levels: 80) 
##                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)            0.38      0.04     0.31     0.45 1.00     4207
## sd(trial)                0.01      0.00     0.00     0.01 1.00     3694
## cor(Intercept,trial)    -0.62      0.12    -0.82    -0.35 1.00     9370
##                      Tail_ESS
## sd(Intercept)            6443
## sd(trial)                4265
## cor(Intercept,trial)     7921
## 
## Population-Level Effects: 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept              3.52      0.06     3.41     3.63 1.00     2993     5184
## conditionexpansive    -0.04      0.07    -0.18     0.10 1.00     3431     5201
## trial                  0.01      0.00     0.00     0.01 1.00     8485     8646
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape     7.67      0.28     7.15     8.22 1.00    17805     8817
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Credible Intervals

model.4.posteriors.CI <- 
model.4.posteriors %>% 
  group_by(condition) %>% 
  median_qi(.epred, width = .95)
model.4.posteriors.CI
condition .epred .epred.lower .epred.upper width width.lower width.upper .width .point .interval
constrictive 39.03590 18.37053 58.48386 0.95 0.95 0.95 0.95 median qi
expansive 37.50873 17.60959 55.62679 0.95 0.95 0.95 0.95 median qi 
model.4.posteriors %>% 
ggplot() +
  geom_histogram(aes(x = .epred, fill = condition), alpha = .5, color = NA) +
  geom_point(model.4.posteriors.CI, 
             mapping = aes(x = .epred, y = 0), size = 3) + 
  geom_errorbarh(model.4.posteriors.CI, 
                 mapping = aes(x = .epred, xmin = .epred.lower, xmax = .epred.upper, y = 0), height = 0, linewidth = 1.5) + 
  facet_wrap(condition ~ ., ncol = 1) + 
  scale_x_continuous(limits = c(0, 80)) + 
  scale_y_continuous(expand = c(.02,.02)) + 
  scale_color_theme() +
  theme_density_x

subtraction

model.4.posteriors_diff <- 
model.4.posteriors %>% 
  ungroup() %>% 
  compare_levels(variable = .epred, by = condition) %>% 
  ungroup()
head(model.4.posteriors_diff)
.draw condition .epred
1 expansive - constrictive 0.4026191
2 expansive - constrictive 0.9758489
3 expansive - constrictive -9.0252287
4 expansive - constrictive -5.1882260
5 expansive - constrictive -2.6665306
6 expansive - constrictive -3.4653370 
model.4.posteriors_diff.CI <-
model.4.posteriors_diff %>% 
  median_qi(.epred)

model.4.posteriors_diff.CI
.epred .lower .upper .width .point .interval
-1.356592 -7.313397 4.011888 0.95 median qi 
model.4.posteriors_diff %>% 
ggplot() +
  geom_histogram(aes(x = .epred), alpha = .5,  fill = theme_blue, color = NA, adjust = 2) +
  geom_point(model.4.posteriors_diff.CI, 
             mapping = aes(x = .epred, y = 0), size = 3) + 
  geom_errorbarh(model.4.posteriors_diff.CI, 
                 mapping = aes(x = .epred, xmin = .lower, xmax = .upper, y = 0), height = 0, linewidth = 1.5) + 
  #scale_x_continuous(limits = c(-50, 50)) + 
  xlab('expansive - constrictive') +
  ggtitle('Mean difference in expansive and constrictive') + 
  geom_vline(xintercept = 0, linetype = 2) + 
  scale_y_continuous(expand = c(.02,.02)) + 
  scale_x_continuous(limits = c(-10, 10), breaks = seq(-10, 10, by = 5)) + 
  scale_color_theme() +
  theme_density_x

Session info

sessionInfo()
## R version 4.2.2 (2022-10-31)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Ventura 13.2
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] cmdstanr_0.5.2      knitr_1.39          ggplot2_3.4.0      
##  [4] bayesplot_1.9.0     ggdist_3.1.1        tidybayes_3.0.2    
##  [7] brms_2.17.0         Rcpp_1.0.8.3        modelr_0.1.8       
## [10] broom.mixed_0.2.9.4 broom_1.0.0         patchwork_1.1.2    
## [13] gtools_3.9.2.1      forcats_0.5.1       tidyr_1.2.0        
## [16] purrr_0.3.4         tibble_3.1.7        dplyr_1.0.9        
## 
## loaded via a namespace (and not attached):
##   [1] colorspace_2.0-3     ellipsis_0.3.2       ggridges_0.5.3      
##   [4] markdown_1.1         base64enc_0.1-3      rstudioapi_0.13     
##   [7] listenv_0.8.0        furrr_0.3.1          farver_2.1.1        
##  [10] rstan_2.21.5         bit64_4.0.5          svUnit_1.0.6        
##  [13] DT_0.23              fansi_1.0.3          mvtnorm_1.1-3       
##  [16] diffobj_0.3.5        bridgesampling_1.1-2 codetools_0.2-18    
##  [19] splines_4.2.2        shinythemes_1.2.0    jsonlite_1.8.0      
##  [22] shiny_1.7.1          readr_2.1.2          compiler_4.2.2      
##  [25] backports_1.4.1      assertthat_0.2.1     Matrix_1.5-1        
##  [28] fastmap_1.1.0        cli_3.4.1            later_1.3.0         
##  [31] htmltools_0.5.2      prettyunits_1.1.1    tools_4.2.2         
##  [34] igraph_1.3.1         coda_0.19-4          gtable_0.3.0        
##  [37] glue_1.6.2           reshape2_1.4.4       posterior_1.2.1     
##  [40] jquerylib_0.1.4      vctrs_0.5.1          nlme_3.1-160        
##  [43] crosstalk_1.2.0      tensorA_0.36.2       xfun_0.31           
##  [46] stringr_1.4.0        globals_0.15.0       ps_1.7.0            
##  [49] mime_0.12            miniUI_0.1.1.1       lifecycle_1.0.3     
##  [52] future_1.26.1        zoo_1.8-10           scales_1.2.0        
##  [55] vroom_1.5.7          colourpicker_1.1.1   hms_1.1.1           
##  [58] promises_1.2.0.1     Brobdingnag_1.2-7    parallel_4.2.2      
##  [61] inline_0.3.19        shinystan_2.6.0      yaml_2.3.5          
##  [64] gridExtra_2.3        loo_2.5.1            StanHeaders_2.21.0-7
##  [67] sass_0.4.1           stringi_1.7.6        highr_0.9           
##  [70] dygraphs_1.1.1.6     checkmate_2.1.0      pkgbuild_1.3.1      
##  [73] rlang_1.0.6          pkgconfig_2.0.3      matrixStats_0.62.0  
##  [76] distributional_0.3.0 evaluate_0.15        lattice_0.20-45     
##  [79] labeling_0.4.2       rstantools_2.2.0     htmlwidgets_1.5.4   
##  [82] cowplot_1.1.1        bit_4.0.4            tidyselect_1.1.2    
##  [85] processx_3.5.3       parallelly_1.31.1    plyr_1.8.7          
##  [88] magrittr_2.0.3       R6_2.5.1             generics_0.1.3      
##  [91] DBI_1.1.2            pillar_1.7.0         withr_2.5.0         
##  [94] xts_0.12.1           abind_1.4-5          crayon_1.5.1        
##  [97] arrayhelpers_1.1-0   utf8_1.2.2           tzdb_0.3.0          
## [100] rmarkdown_2.14       grid_4.2.2           data.table_1.14.2   
## [103] callr_3.7.0          threejs_0.3.3        digest_0.6.29       
## [106] xtable_1.8-4         httpuv_1.6.5         RcppParallel_5.1.5  
## [109] stats4_4.2.2         munsell_0.5.0        bslib_0.3.1         
## [112] shinyjs_2.1.0